# -*- coding: utf-8 -*-
###################################################
#    This file is part of blockIMH.
#
#    blockIMH is free software: you can redistribute it and/or modify
#    it under the terms of the GNU General Public License as published by
#    the Free Software Foundation, either version 3 of the License, or
#    (at your option) any later version.
#
#    blockIMH is distributed in the hope that it will be useful,
#    but WITHOUT ANY WARRANTY; without even the implied warranty of
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#    GNU General Public License for more details.
#
#    You should have received a copy of the GNU General Public License
#    along with blockIMH.  If not, see <http://www.gnu.org/licenses/>.
###################################################
#! /usr/bin/env python

from __future__ import division
import os
import sys
import scipy.weave as weave
from scipy.stats import norm, cauchy, uniform
from scipy import special
from numpy import random, exp, sqrt, power, log, newaxis, zeros, \
        zeros_like, float32, int32, repeat, minimum, ones, mean, \
        average, sum, var, arange
try:
    import pycuda
    import pycuda.autoinit
    import pycuda.gpuarray as gpuarray
    import pycuda.driver as cuda
    import pycuda.cumath as cumath
    from pycuda.compiler import SourceModule
    from mycudafunctions import *
except:
    CUDA_ENABLED = False
from IMHprobit import blockIMH
import probit
import permutations
from localfolder import get_path
import rpy2.robjects as robjects

def h(x):
    return (x)

## Nb of parallel processors
try:
    for arg in sys.argv:
        arg = arg.replace("\s+", "")
        if "p=" in arg:
            p = int(arg.replace("--p=", ""))
        if "b=" in arg:
            b = int(arg.replace("--b=", ""))
        if "nbsim=" in arg:
            nbsim = int(arg.replace("--nbsim=", ""))
except:
    raise NameError("""
problem processing the arguments: the only available options are:
--p=X where X is an integer denoting a block's length,
--b=X where X is an integer denoting the number of blocks (hence T = b * p is the chain's length),
--nbsim=X where X is an integer denoting the number of independent runs to perform.
(don't put spaces around the "=" signs)
""")
if not vars().has_key("p"):
    print "no p given in argument, default is p = 4"
    p = 4
if not vars().has_key("b"):
    print "no b given in argument, default is b = 1"
    b = 1
if not vars().has_key("nbsim"):
    print "no nbsim given in argument, default is nbsim = 100"
    nbsim = 100

T = p * b
print """Starting block IMH algorithm on the probit example with:
    p = %i
    b = %i
    T = %i
    %i independent runs
""" % (p, b, T, nbsim)
x_0 = probit.estimMLE
print "starting point:"
print x_0
basefolder = get_path()
resultsfolder = os.path.join(basefolder, "results")
barplotname = os.path.join(resultsfolder, "probitbarplotN%iP%iT%ic%i.eps" % (nbsim, p, T, probit.scalefactor))
Rfilename = os.path.join(resultsfolder, "probitresultsN%iP%iT%ic%i.R" % (nbsim, p, T, probit.scalefactor))
#barplotname = os.path.join(resultsfolder, "probitTemp.eps")
#Rfilename = os.path.join(resultsfolder, "probitTemp.R")

parameterdimension = x_0.shape[0]

results = {}
results["basicestimvector"] = zeros((nbsim, parameterdimension))
results["blockestimvectorRandPermut"] = zeros((nbsim, parameterdimension))
results["estimRB1vectorRandPermut"] = zeros((nbsim, parameterdimension))
results["estimRB2vectorRandPermut"] = zeros((nbsim, parameterdimension))
results["estimIS"] = zeros((nbsim, parameterdimension))
#permutationmatrix = permutations.createRandomPermutation(p)

#imh = blockIMH(x_0, T, p, h, permutationmatrix)
##print imh
#print "results"
#print "basic", imh["basicestim"]
#print "block", imh["blockestimRandPermut"]
#print "RB1  ", imh["estimRB1RandPermut"]
#print "RB2  ", imh["estimRB2RandPermut"]
#indices = arange(start = 0, stop = T, step = 1)
#robjects.r(""" par(mfrow = c(%i, 2)) """ % parameterdimension)
#for i in range(parameterdimension):
#    x = imh["chain"][indices, i]
#    rx = robjects.FloatVector(x)
#    robjects.r.plot(rx, xlab = "", ylab = "", type ="l")
#    robjects.r.hist(rx, xlab = "", ylab = "", main = "", col = "red", prob = 1, nclass = 100)
#raw_input()

import cProfile
cProfile.run("""
for simu in range(nbsim):
    print simu
    imh = blockIMH(x_0, T, p, h)
    results["basicestimvector"][simu, :] = imh["basicestim"]
    results["blockestimvectorRandPermut"][simu, :] = imh["blockestimRandPermut"]
    results["estimRB1vectorRandPermut"][simu, :] = imh["estimRB1RandPermut"]
    results["estimRB2vectorRandPermut"][simu, :] = imh["estimRB2RandPermut"]
    results["estimIS"][simu] = imh["estimIS"]
""", "prof")

import pstats
pS = pstats.Stats('prof')
pS.sort_stats("cumulative").print_stats(10)
pS.sort_stats("time").print_stats(10)
pS.sort_stats("cumulative").print_stats("omega")

import cPickle
picklefilename = "results/results-probit-p%iT%in%i-scale%i.cpickle" % (p, T, nbsim, probit.scalefactor)
f = open(picklefilename, "w")
cPickle.dump(results, f)
f.close()
import cPickle
import pickletoRdata
pickletoRdata.pickle2RData(picklefilename)
os.remove(picklefilename)

print "Goodbye!"




